{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {},
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D1_BayesianDecisions/student/W3D1_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a> &nbsp; <a href=\"https://kaggle.com/kernels/welcome?src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W7_BayesianDecisions/student/W7_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://kaggle.com/static/images/open-in-kaggle.svg\" alt=\"Open in Kaggle\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Tutorial 2: Bayesian inference and decisions with continuous hidden state\n",
    "\n",
    "__Content creators:__ Eric DeWitt, Xaq Pitkow, Saeed Salehi, Ella Batty\n",
    "\n",
    "__Content modified:__ Kai Chen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Tutorial Objectives\n",
    "\n",
    "This notebook introduces you to Gaussians and Bayes' rule for continuous distributions, allowing us to model simple put powerful combinations of prior information and new measurements. In this notebook you will work through the same ideas we explored in the binary state tutorial, but you will be introduced to a new problem: finding and guiding Astrocat! You will see this problem again in more complex ways in the following days.\n",
    "\n",
    "In this notebook, you will:\n",
    "\n",
    "1. Learn about the Gaussian distribution and its nice properies\n",
    "2. Explore how we can extend the ideas from the binary hidden tutorial to continuous distributions\n",
    "3. Explore how different priors can produce more complex posteriors.\n",
    "4. Explore Loss functions often used in inference and complex utility functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Setup  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import multivariate_normal\n",
    "from scipy.stats import gamma as gamma_distribution\n",
    "from matplotlib.transforms import Affine2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output\n",
    "from ipywidgets import FloatSlider, Dropdown, interactive_output\n",
    "from ipywidgets import interact, fixed, HBox, Layout, VBox, interactive, Label\n",
    "\n",
    "nma_style = {\n",
    "    'figure.figsize' : (8, 6),\n",
    "    'figure.autolayout' : True,\n",
    "    'font.size' : 15,\n",
    "    'xtick.labelsize' : 'small',\n",
    "    'ytick.labelsize' : 'small',\n",
    "    'legend.fontsize' : 'small',\n",
    "    'axes.spines.top' : False,\n",
    "    'axes.spines.right' : False,\n",
    "    'xtick.major.size' : 5,\n",
    "    'ytick.major.size' : 5,\n",
    "}\n",
    "for key, value in nma_style.items():\n",
    "    plt.rcParams[key] = value\n",
    "\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title Plotting Functions\n",
    "\n",
    "def plot_mixture_prior(x, gaussian1, gaussian2, combined):\n",
    "    \"\"\" Plots a prior made of a mixture of gaussians\n",
    "\n",
    "    Args:\n",
    "      x (numpy array of floats):         points at which the likelihood has been evaluated\n",
    "      gaussian1 (numpy array of floats): normalized probabilities for Gaussian 1 evaluated at each `x`\n",
    "      gaussian2 (numpy array of floats): normalized probabilities for Gaussian 2 evaluated at each `x`\n",
    "      posterior (numpy array of floats): normalized probabilities for the posterior evaluated at each `x`\n",
    "\n",
    "    Returns:\n",
    "      Nothing\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(x, gaussian1, '--b', LineWidth=2, label='Gaussian 1')\n",
    "    ax.plot(x, gaussian2, '-.b', LineWidth=2, label='Gaussian 2')\n",
    "    ax.plot(x, combined, '-r', LineWidth=2, label='Gaussian Mixture')\n",
    "    ax.legend()\n",
    "    ax.set_ylabel('Probability')\n",
    "    ax.set_xlabel('Orientation (Degrees)')\n",
    "\n",
    "\n",
    "def plot_gaussian(μ, σ):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    y = gaussian(x, μ, σ)\n",
    "\n",
    "    plt.figure(figsize=(6, 4))\n",
    "    plt.plot(x, y, c='blue')\n",
    "    plt.fill_between(x, y, color='b', alpha=0.2)\n",
    "    plt.ylabel('$\\mathcal{N}(x, \\mu, \\sigma^2)$')\n",
    "    plt.xlabel('x')\n",
    "    plt.yticks([])\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_losses(μ, σ):\n",
    "    x = np.linspace(-2, 2, 400, endpoint=True)\n",
    "    y = gaussian(x, μ, σ)\n",
    "    error = x - μ\n",
    "\n",
    "    mse_loss = (error)**2\n",
    "    abs_loss = np.abs(error)\n",
    "    zero_one_loss = (np.abs(error) >= 0.02).astype(np.float)\n",
    "\n",
    "    fig, (ax_gaus, ax_error) = plt.subplots(2, 1, figsize=(6, 8))\n",
    "    ax_gaus.plot(x, y, color='blue', label='true distribution')\n",
    "    ax_gaus.fill_between(x, y, color='blue', alpha=0.2)\n",
    "    ax_gaus.set_ylabel('$\\\\mathcal{N}(x, \\\\mu, \\\\sigma^2)$')\n",
    "    ax_gaus.set_xlabel('x')\n",
    "    ax_gaus.set_yticks([])\n",
    "    ax_gaus.legend(loc='upper right')\n",
    "\n",
    "    ax_error.plot(x, mse_loss, color='c', label='Mean Squared Error', linewidth=3)\n",
    "    ax_error.plot(x, abs_loss, color='m', label='Absolute Error', linewidth=3)\n",
    "    ax_error.plot(x, zero_one_loss, color='y', label='Zero-One Loss', linewidth=3)\n",
    "    ax_error.legend(loc='upper right')\n",
    "    ax_error.set_xlabel('$\\\\hat{\\\\mu}$')\n",
    "    ax_error.set_ylabel('Error')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def gaussian_mixture(mu1, mu2, sigma1, sigma2, factor):\n",
    "    assert 0.0 < factor < 1.0\n",
    "    x = np.linspace(-7.0, 7.0, 1000, endpoint=True)\n",
    "    y_1 = gaussian(x, mu1, sigma1)\n",
    "    y_2 = gaussian(x, mu2, sigma2)\n",
    "    mixture = y_1 * factor + y_2 * (1.0 - factor)\n",
    "\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.plot(x, y_1, c='deepskyblue', label='p(x)', linewidth=3.0)\n",
    "    plt.fill_between(x, y_1, color='deepskyblue', alpha=0.2)\n",
    "    plt.plot(x, y_2, c='aquamarine', label='q(x)', linewidth=3.0)\n",
    "    plt.fill_between(x, y_2, color='aquamarine', alpha=0.2)\n",
    "    plt.plot(x, mixture, c='b', label='$\\pi \\cdot p(x) + (1-\\pi) \\cdot q(x)$',  linewidth=3.0)\n",
    "    plt.fill_between(x, mixture, color='b', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    # plt.ylabel('$f(x)$')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_utility_mixture_dist(mu1, sigma1, mu2, sigma2, mu_g, sigma_g,\n",
    "                              mu_loc, mu_dist, plot_utility_row=True):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "    mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    posterior = gaussian(x, mu_post, sigma_post)\n",
    "    gain = gaussian(x, mu_g, sigma_g)/2\n",
    "\n",
    "    sigma_mix, factor = 1.0, 0.5\n",
    "    mu_mix1, mu_mix2 = mu_loc - mu_dist/2, mu_loc + mu_dist/2\n",
    "    gaus_mix1, gaus_mix2 = gaussian(x, mu_mix1, sigma_mix), gaussian(x, mu_mix2, sigma_mix)\n",
    "    loss = factor * gaus_mix1 + (1 - factor) * gaus_mix2\n",
    "    utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "\n",
    "    if plot_utility_row:\n",
    "        plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)\n",
    "    else:\n",
    "        plot_bayes_row(x, prior, likelihood, posterior)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_mvn2d(mu1, mu2, sigma1, sigma2, corr):\n",
    "    x, y = np.mgrid[-2:2:.02, -2:2:.02]\n",
    "    cov12 = corr * sigma1 * sigma2\n",
    "    z = mvn2d(x, y, mu1, mu2, sigma1, sigma2, cov12)\n",
    "\n",
    "    plt.figure(figsize=(6, 6))\n",
    "    plt.contourf(x, y, z, cmap='Reds')\n",
    "    plt.axis(\"off\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_marginal(sigma1, sigma2, c_x, c_y, corr):\n",
    "    mu1, mu2 = 0.0, 0.0\n",
    "    cov12 = corr * sigma1 * sigma2\n",
    "    xx, yy = np.mgrid[-2:2:.02, -2:2:.02]\n",
    "    x, y = xx[:, 0], yy[0]\n",
    "    p_x = gaussian(x, mu1, sigma1)\n",
    "    p_y = gaussian(y, mu2, sigma2)\n",
    "    zz = mvn2d(xx, yy, mu1, mu2, sigma1, sigma2, cov12)\n",
    "\n",
    "    mu_x_y = mu1+cov12*(c_y-mu2)/sigma2**2\n",
    "    mu_y_x = mu2+cov12*(c_x-mu1)/sigma1**2\n",
    "    sigma_x_y = np.sqrt(sigma2**2 - cov12**2/sigma1**2)\n",
    "    sigma_y_x = np.sqrt(sigma1**2-cov12**2/sigma2**2)\n",
    "    p_x_y = gaussian(x, mu_x_y, sigma_x_y)\n",
    "    p_y_x = gaussian(x, mu_y_x, sigma_y_x)\n",
    "\n",
    "    p_c_y = gaussian(mu_x_y-sigma_x_y, mu_x_y, sigma_x_y)\n",
    "    p_c_x = gaussian(mu_y_x-sigma_y_x, mu_y_x, sigma_y_x)\n",
    "\n",
    "    # definitions for the axes\n",
    "    left, width = 0.1, 0.65\n",
    "    bottom, height = 0.1, 0.65\n",
    "    spacing = 0.01\n",
    "\n",
    "    rect_z = [left, bottom, width, height]\n",
    "    rect_x = [left, bottom + height + spacing, width, 0.2]\n",
    "    rect_y = [left + width + spacing, bottom, 0.2, height]\n",
    "\n",
    "    # start with a square Figure\n",
    "    fig = plt.figure(figsize=(8, 8))\n",
    "\n",
    "    ax_z = fig.add_axes(rect_z)\n",
    "    ax_x = fig.add_axes(rect_x, sharex=ax_z)\n",
    "    ax_y = fig.add_axes(rect_y, sharey=ax_z)\n",
    "\n",
    "    ax_z.set_axis_off()\n",
    "    ax_x.set_axis_off()\n",
    "    ax_y.set_axis_off()\n",
    "    ax_x.set_xlim(np.min(x), np.max(x))\n",
    "    ax_y.set_ylim(np.min(y), np.max(y))\n",
    "\n",
    "    ax_z.contourf(xx, yy, zz, cmap='Greys')\n",
    "    ax_z.hlines(c_y, mu_x_y-sigma_x_y, mu_x_y+sigma_x_y, color='c', zorder=9, linewidth=3)\n",
    "    ax_z.vlines(c_x, mu_y_x-sigma_y_x, mu_y_x+sigma_y_x, color='m', zorder=9, linewidth=3)\n",
    "\n",
    "    ax_x.plot(x, p_x, label='$p(x)$', c = 'b', linewidth=3)\n",
    "    ax_x.plot(x, p_x_y, label='$p(x|y = C_y)$', c = 'c', linestyle='dashed', linewidth=3)\n",
    "    ax_x.hlines(p_c_y, mu_x_y-sigma_x_y, mu_x_y+sigma_x_y, color='c', linestyle='dashed', linewidth=3)\n",
    "\n",
    "    ax_y.plot(p_y, y, label='$p(y)$', c = 'r', linewidth=3)\n",
    "    ax_y.plot(p_y_x, y, label='$p(y|x = C_x)$', c = 'm', linestyle='dashed', linewidth=3)\n",
    "    ax_y.vlines(p_c_x, mu_y_x-sigma_y_x, mu_y_x+sigma_y_x, color='m', linestyle='dashed', linewidth=3)\n",
    "\n",
    "    ax_x.legend(loc=\"upper left\", frameon=False)\n",
    "    ax_y.legend(loc=\"lower right\", frameon=False)\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_bayes(mu1, mu2, sigma1, sigma2):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "\n",
    "    mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    posterior = gaussian(x, mu_post, sigma_post)\n",
    "\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.plot(x, prior, c='b', label='prior')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='likelihood')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    plt.plot(x, posterior, c='k', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.ylabel('$\\mathcal{N}(x, \\mu, \\sigma^2)$')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()\n",
    "\n",
    "def plot_information(mu1, sigma1, mu2, sigma2):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    mu3, sigma3 = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "    posterior = gaussian(x, mu3, sigma3)\n",
    "\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.plot(x, prior, c='b', label='Satellite')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='Space Mouse')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    plt.plot(x, posterior, c='k', label='Center')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.ylabel('$\\mathcal{N}(x, \\mu, \\sigma^2)$')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_information_global(mu3, sigma3, mu1, mu2):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    sigma1, sigma2 = reverse_product(mu3, sigma3, mu1, mu2)\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "    posterior = gaussian(x, mu3, sigma3)\n",
    "\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.plot(x, prior, c='b', label='Satellite')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='Space Mouse')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    plt.plot(x, posterior, c='k', label='Center')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.ylabel('$\\mathcal{N}(x, \\mu, \\sigma^2)$')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_loss_utility_gaussian(loss_f, mu, sigma, mu_true):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    posterior = gaussian(x, mu, sigma)\n",
    "    y_label = \"$p(x)$\"\n",
    "\n",
    "    plot_loss_utility(x, posterior, loss_f, mu_true, y_label)\n",
    "\n",
    "\n",
    "def plot_loss_utility_mixture(loss_f, mu1, mu2, sigma1, sigma2, factor, mu_true):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    y_1 = gaussian(x, mu1, sigma1)\n",
    "    y_2 = gaussian(x, mu2, sigma2)\n",
    "    posterior = y_1 * factor + y_2 * (1.0 - factor)\n",
    "    y_label = \"$\\pi \\cdot p(x) + (1-\\pi) \\cdot q(x)$\"\n",
    "    plot_loss_utility(x, posterior, loss_f, mu_true, y_label)\n",
    "\n",
    "\n",
    "def plot_loss_utility(x, posterior, loss_f, mu_true, y_label):\n",
    "    mean, median, mode = calc_mean_mode_median(x, posterior)\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    utility = calc_expected_loss(loss_f, posterior, x)\n",
    "    min_expected_loss = x[np.argmin(utility)]\n",
    "\n",
    "    plt.figure(figsize=(12, 8))\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.title(\"Probability\")\n",
    "    plt.plot(x, posterior, c='b')\n",
    "    plt.fill_between(x, posterior, color='b', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel(y_label)\n",
    "    plt.axvline(mean, ls='dashed', color='red', label='Mean')\n",
    "    plt.axvline(median, ls='dashdot', color='blue', label='Median')\n",
    "    plt.axvline(mode, ls='dotted', color='green', label='Mode')\n",
    "    plt.legend(loc=\"upper left\")\n",
    "\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.title(loss_f)\n",
    "    plt.plot(x, loss, c='c', label=loss_f)\n",
    "    # plt.fill_between(x, loss, color='c', alpha=0.2)\n",
    "    plt.ylabel('loss')\n",
    "    # plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.title(\"Expected Loss\")\n",
    "    plt.plot(x, utility, c='y', label='$\\mathbb{E}[L]$')\n",
    "    plt.axvline(min_expected_loss, ls='dashed', color='red', label='$Min~ \\mathbb{E}[Loss]$')\n",
    "    # plt.fill_between(x, utility, color='y', alpha=0.2)\n",
    "    plt.legend(loc=\"lower right\")\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('$\\mathbb{E}[L]$')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_loss_utility_bayes(mu1, mu2, sigma1, sigma2, mu_true, loss_f):\n",
    "    x = np.linspace(-4, 4, 1000, endpoint=True)\n",
    "\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "\n",
    "    mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    posterior = gaussian(x, mu_post, sigma_post)\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    utility = - calc_expected_loss(loss_f, posterior, x)\n",
    "\n",
    "    plt.figure(figsize=(18, 5))\n",
    "    plt.subplot(1, 3, 1)\n",
    "\n",
    "    plt.title(\"Posterior distribution\")\n",
    "    plt.plot(x, prior, c='b', label='prior')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='likelihood')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    plt.plot(x, posterior, c='k', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    # plt.ylabel('$f(x)$')\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(1, 3, 2)\n",
    "    plt.title(loss_f)\n",
    "    plt.plot(x, loss, c='c')\n",
    "    # plt.fill_between(x, loss, color='c', alpha=0.2)\n",
    "    plt.ylabel('loss')\n",
    "\n",
    "    plt.subplot(1, 3, 3)\n",
    "    plt.title(\"Expected utility\")\n",
    "    plt.plot(x, utility, c='y', label='utility')\n",
    "    # plt.fill_between(x, utility, color='y', alpha=0.2)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_simple_utility_gaussian(mu, sigma, mu_g, mu_c, sigma_g, sigma_c):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    posterior = gaussian(x, mu, sigma)\n",
    "    gain = gaussian(x, mu_g, sigma_g)\n",
    "    loss = gaussian(x, mu_c, sigma_c)\n",
    "    utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "\n",
    "    plt.figure(figsize=(15, 5))\n",
    "    plt.subplot(1, 3, 1)\n",
    "    plt.title(\"Probability\")\n",
    "    plt.plot(x, posterior, c='b', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='b', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    # plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(1, 3, 2)\n",
    "    plt.title(\"utility function\")\n",
    "    plt.plot(x, gain, c='m', label='gain')\n",
    "    # plt.fill_between(x, gain, color='m', alpha=0.2)\n",
    "    plt.plot(x, -loss, c='c', label='loss')\n",
    "    # plt.fill_between(x, -loss, color='c', alpha=0.2)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "\n",
    "    plt.subplot(1, 3, 3)\n",
    "    plt.title(\"expected utility\")\n",
    "    plt.plot(x, utility, c='y', label='utility')\n",
    "    # plt.fill_between(x, utility, color='y', alpha=0.2)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_utility_gaussian(mu1, mu2, sigma1, sigma2, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "\n",
    "    mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    posterior = gaussian(x, mu_post, sigma_post)\n",
    "\n",
    "    if plot_utility_row:\n",
    "        gain = gaussian(x, mu_g, sigma_g)\n",
    "        loss = gaussian(x, mu_c, sigma_c)\n",
    "        utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "        plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)\n",
    "    else:\n",
    "        plot_bayes_row(x, prior, likelihood, posterior)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_utility_mixture(mu_m1, mu_m2, sigma_m1, sigma_m2, factor,\n",
    "                         mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    y_1 = gaussian(x, mu_m1, sigma_m1)\n",
    "    y_2 = gaussian(x, mu_m2, sigma_m2)\n",
    "    prior = y_1 * factor + y_2 * (1.0 - factor)\n",
    "\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "\n",
    "    posterior = np.multiply(prior, likelihood)\n",
    "    posterior = posterior / (posterior.sum() * (x[1] - x[0]))\n",
    "\n",
    "    if plot_utility_row:\n",
    "        gain = gaussian(x, mu_g, sigma_g)\n",
    "        loss = gaussian(x, mu_c, sigma_c)\n",
    "        utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "        plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)\n",
    "    else:\n",
    "        plot_bayes_row(x, prior, likelihood, posterior)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_utility_uniform(mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    prior = np.ones_like(x) / (x.max() - x.min())\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "\n",
    "    posterior = likelihood\n",
    "    # posterior = np.multiply(prior, likelihood)\n",
    "    # posterior = posterior / (posterior.sum() * (x[1] - x[0]))\n",
    "\n",
    "    if plot_utility_row:\n",
    "        gain = gaussian(x, mu_g, sigma_g)\n",
    "        loss = gaussian(x, mu_c, sigma_c)\n",
    "        utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "        plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)\n",
    "    else:\n",
    "        plot_bayes_row(x, prior, likelihood, posterior)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_utility_gamma(alpha, beta, offset, mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):\n",
    "    x = np.linspace(-4, 10, 1000, endpoint=True)\n",
    "    prior = gamma_pdf(x-offset, alpha, beta)\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "\n",
    "    posterior = np.multiply(prior, likelihood)\n",
    "    posterior = posterior / (posterior.sum() * (x[1] - x[0]))\n",
    "\n",
    "    if plot_utility_row:\n",
    "        gain = gaussian(x, mu_g, sigma_g)\n",
    "        loss = gaussian(x, mu_c, sigma_c)\n",
    "        utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)\n",
    "        plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)\n",
    "    else:\n",
    "        plot_bayes_row(x, prior, likelihood, posterior)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_bayes_row(x, prior, likelihood, posterior):\n",
    "\n",
    "    mean, median, mode = calc_mean_mode_median(x, posterior)\n",
    "\n",
    "    plt.figure(figsize=(12, 4))\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.title(\"Prior and likelihood distribution\")\n",
    "    plt.plot(x, prior, c='b', label='prior')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='likelihood')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    # plt.plot(x, posterior, c='k', label='posterior')\n",
    "    # plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    # plt.ylabel('$f(x)$')\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.title(\"Posterior distribution\")\n",
    "    plt.plot(x, posterior, c='k', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.1)\n",
    "    plt.axvline(mean, ls='dashed', color='red', label='Mean')\n",
    "    plt.axvline(median, ls='dashdot', color='blue', label='Median')\n",
    "    plt.axvline(mode, ls='dotted', color='green', label='Mode')\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.yticks([])\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility):\n",
    "\n",
    "    mean, median, mode = calc_mean_mode_median(x, posterior)\n",
    "    max_utility = x[np.argmax(utility)]\n",
    "\n",
    "    plt.figure(figsize=(12, 8))\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.title(\"Prior and likelihood distribution\")\n",
    "    plt.plot(x, prior, c='b', label='prior')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='likelihood')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    # plt.plot(x, posterior, c='k', label='posterior')\n",
    "    # plt.fill_between(x, posterior, color='k', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    # plt.ylabel('$f(x)$')\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.title(\"Posterior distribution\")\n",
    "    plt.plot(x, posterior, c='k', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.1)\n",
    "    plt.axvline(mean, ls='dashed', color='red', label='Mean')\n",
    "    plt.axvline(median, ls='dashdot', color='blue', label='Median')\n",
    "    plt.axvline(mode, ls='dotted', color='green', label='Mode')\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.yticks([])\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.title(\"utility function\")\n",
    "    plt.plot(x, gain, c='m', label='gain')\n",
    "    # plt.fill_between(x, gain, color='m', alpha=0.2)\n",
    "    plt.plot(x, -loss, c='c', label='loss')\n",
    "    # plt.fill_between(x, -loss, color='c', alpha=0.2)\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.title(\"expected utility\")\n",
    "    plt.plot(x, utility, c='y', label='utility')\n",
    "    # plt.fill_between(x, utility, color='y', alpha=0.2)\n",
    "    plt.axvline(max_utility, ls='dashed', color='red', label='Max utility')\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('utility')\n",
    "    plt.legend(loc=\"lower right\")\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_bayes_loss_utility_gaussian(loss_f, mu_true, mu1, mu2, sigma1, sigma2):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "\n",
    "    prior = gaussian(x, mu1, sigma1)\n",
    "    likelihood = gaussian(x, mu2, sigma2)\n",
    "    mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)\n",
    "    posterior = gaussian(x, mu_post, sigma_post)\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_bayes_loss_utility_uniform(loss_f, mu_true, mu, sigma):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "\n",
    "    prior = np.ones_like(x) / (x.max() - x.min())\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "    posterior = likelihood\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_bayes_loss_utility_gamma(loss_f, mu_true, alpha, beta, offset, mu, sigma):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    prior = gamma_pdf(x-offset, alpha, beta)\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "    posterior = np.multiply(prior, likelihood)\n",
    "    posterior = posterior / (posterior.sum() * (x[1] - x[0]))\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_bayes_loss_utility_mixture(loss_f, mu_true, mu_m1, mu_m2, sigma_m1, sigma_m2, factor, mu, sigma):\n",
    "    x = np.linspace(-7, 7, 1000, endpoint=True)\n",
    "    y_1 = gaussian(x, mu_m1, sigma_m1)\n",
    "    y_2 = gaussian(x, mu_m2, sigma_m2)\n",
    "    prior = y_1 * factor + y_2 * (1.0 - factor)\n",
    "    likelihood = gaussian(x, mu, sigma)\n",
    "\n",
    "    posterior = np.multiply(prior, likelihood)\n",
    "    posterior = posterior / (posterior.sum() * (x[1] - x[0]))\n",
    "\n",
    "    loss = calc_loss_func(loss_f, mu_true, x)\n",
    "\n",
    "    plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "def plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f):\n",
    "\n",
    "    mean, median, mode = calc_mean_mode_median(x, posterior)\n",
    "    expected_loss = calc_expected_loss(loss_f, posterior, x)\n",
    "    min_expected_loss = x[np.argmin(expected_loss)]\n",
    "\n",
    "    plt.figure(figsize=(12, 8))\n",
    "\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.title(\"Prior and Likelihood\")\n",
    "    plt.plot(x, prior, c='b', label='prior')\n",
    "    plt.fill_between(x, prior, color='b', alpha=0.2)\n",
    "    plt.plot(x, likelihood, c='r', label='likelihood')\n",
    "    plt.fill_between(x, likelihood, color='r', alpha=0.2)\n",
    "    plt.yticks([])\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.title(\"Posterior\")\n",
    "    plt.plot(x, posterior, c='k', label='posterior')\n",
    "    plt.fill_between(x, posterior, color='k', alpha=0.1)\n",
    "    plt.axvline(mean, ls='dashed', color='red', label='Mean')\n",
    "    plt.axvline(median, ls='dashdot', color='blue', label='Median')\n",
    "    plt.axvline(mode, ls='dotted', color='green', label='Mode')\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.yticks([])\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.title(loss_f)\n",
    "    plt.plot(x, loss, c='c', label=loss_f)\n",
    "    # plt.fill_between(x, loss, color='c', alpha=0.2)\n",
    "    plt.ylabel('loss')\n",
    "    plt.xlabel('x')\n",
    "\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.title(\"expected loss\")\n",
    "    plt.plot(x, expected_loss, c='y', label='$\\mathbb{E}[L]$')\n",
    "    # plt.fill_between(x, expected_loss, color='y', alpha=0.2)\n",
    "    plt.axvline(min_expected_loss, ls='dashed', color='red', label='$Min~ \\mathbb{E}[Loss]$')\n",
    "    plt.legend(loc=\"lower right\")\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('$\\mathbb{E}[L]$')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "global global_loss_plot_switcher\n",
    "global_loss_plot_switcher = False\n",
    "def loss_plot_switcher(what_to_plot):\n",
    "    global global_loss_plot_switcher\n",
    "    if global_loss_plot_switcher:\n",
    "        clear_output()\n",
    "    else:\n",
    "        global_loss_plot_switcher = True\n",
    "    if what_to_plot == \"Gaussian\":\n",
    "        loss_f_options = Dropdown(\n",
    "                    options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "                    value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_estimate\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_estimate\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-3.0, max=3.0, step=0.01, value=0.0, description=\"µ_true\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([loss_f_options, mu_true_slider]), VBox([mu_slider, sigma_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_loss_utility_gaussian,\n",
    "                                        {'loss_f': loss_f_options,\n",
    "                                        'mu': mu_slider,\n",
    "                                        'sigma': sigma_slider,\n",
    "                                        'mu_true': mu_true_slider})\n",
    "\n",
    "    elif what_to_plot == \"Mixture of Gaussians\":\n",
    "        loss_f_options = Dropdown(\n",
    "                options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "                value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "\n",
    "        mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_est_p\", continuous_update=True)\n",
    "        mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_est_q\", continuous_update=True)\n",
    "        sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_est_p\", continuous_update=True)\n",
    "        sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_est_q\", continuous_update=True)\n",
    "        factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description=\"π\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-3.0, max=3.0, step=0.01, value=0.0, description=\"µ_true\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([loss_f_options, factor_slider, mu_true_slider]),\n",
    "                          VBox([mu1_slider, sigma1_slider]),\n",
    "                          VBox([mu2_slider, sigma2_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_loss_utility_mixture,\n",
    "                                        {'mu1': mu1_slider,\n",
    "                                        'mu2': mu2_slider,\n",
    "                                        'sigma1': sigma1_slider,\n",
    "                                        'sigma2': sigma2_slider,\n",
    "                                        'factor': factor_slider,\n",
    "                                        'mu_true': mu_true_slider,\n",
    "                                        'loss_f': loss_f_options})\n",
    "    display(widget_ui, widget_out)\n",
    "\n",
    "\n",
    "global global_plot_prior_switcher\n",
    "global_plot_prior_switcher = False\n",
    "def plot_prior_switcher(what_to_plot):\n",
    "    global global_plot_prior_switcher\n",
    "    if global_plot_prior_switcher:\n",
    "        clear_output()\n",
    "    else:\n",
    "        global_plot_prior_switcher = True\n",
    "\n",
    "    if what_to_plot == \"Gaussian\":\n",
    "        mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_prior\", continuous_update=True)\n",
    "        mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_prior\", continuous_update=True)\n",
    "        sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([mu1_slider, sigma1_slider]),\n",
    "                    VBox([mu2_slider, sigma2_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_utility_gaussian,\n",
    "                                        {'mu1': mu1_slider,\n",
    "                                         'mu2': mu2_slider,\n",
    "                                         'sigma1': sigma1_slider,\n",
    "                                         'sigma2': sigma2_slider,\n",
    "                                         'mu_g': fixed(1.0),\n",
    "                                         'mu_c': fixed(-1.0),\n",
    "                                         'sigma_g': fixed(0.5),\n",
    "                                         'sigma_c': fixed(value=0.5),\n",
    "                                         'plot_utility_row': fixed(False)})\n",
    "\n",
    "    elif what_to_plot == \"Uniform\":\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "\n",
    "        widget_ui = VBox([mu_slider, sigma_slider])\n",
    "\n",
    "        widget_out = interactive_output(plot_utility_uniform,\n",
    "                                        {'mu': mu_slider,\n",
    "                                         'sigma': sigma_slider,\n",
    "                                         'mu_g': fixed(1.0),\n",
    "                                         'mu_c': fixed(-1.0),\n",
    "                                         'sigma_g': fixed(0.5),\n",
    "                                         'sigma_c': fixed(value=0.5),\n",
    "                                         'plot_utility_row': fixed(False)})\n",
    "\n",
    "    elif what_to_plot == \"Gamma\":\n",
    "        alpha_slider = FloatSlider(min=1.0, max=10.0, step=0.1, value=2.0, description=\"α_prior\", continuous_update=True)\n",
    "        beta_slider = FloatSlider(min=0.5, max=2.0, step=0.01, value=1.0, description=\"β_prior\", continuous_update=True)\n",
    "        # offset_slider = FloatSlider(min=-6.0, max=2.0, step=0.1, value=0.0, description=\"offset\", continuous_update=True)\n",
    "        offset_slider = fixed(0.0)\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "        gaus_label = Label(value=\"normal likelihood\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "        gamma_label = Label(value=\"gamma prior\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "        widget_ui = HBox([VBox([gamma_label, alpha_slider, beta_slider]),\n",
    "                          VBox([gaus_label, mu_slider, sigma_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_utility_gamma,\n",
    "                                        {'alpha': alpha_slider,\n",
    "                                         'beta': beta_slider,\n",
    "                                         'offset': offset_slider,\n",
    "                                         'mu': mu_slider,\n",
    "                                         'sigma': sigma_slider,\n",
    "                                         'mu_g': fixed(1.0),\n",
    "                                         'mu_c': fixed(-1.0),\n",
    "                                         'sigma_g': fixed(0.5),\n",
    "                                         'sigma_c': fixed(value=0.5),\n",
    "                                         'plot_utility_row': fixed(False)})\n",
    "\n",
    "    elif what_to_plot == \"Mixture of Gaussians\":\n",
    "        mu_m1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_mix_p\", continuous_update=True)\n",
    "        mu_m2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_mix_q\", continuous_update=True)\n",
    "        sigma_m1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_mix_p\", continuous_update=True)\n",
    "        sigma_m2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_mix_q\", continuous_update=True)\n",
    "        factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description=\"π\", continuous_update=True)\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([mu_m1_slider, sigma_m1_slider, factor_slider]),\n",
    "                          VBox([mu_m2_slider, sigma_m2_slider]),\n",
    "                          VBox([mu_slider, sigma_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_utility_mixture,\n",
    "                                        {'mu_m1': mu_m1_slider,\n",
    "                                         'mu_m2': mu_m2_slider,\n",
    "                                         'sigma_m1': sigma_m1_slider,\n",
    "                                         'sigma_m2': sigma_m2_slider,\n",
    "                                         'factor': factor_slider,\n",
    "                                         'mu': mu_slider,\n",
    "                                         'sigma': sigma_slider,\n",
    "                                         'mu_g': fixed(1.0),\n",
    "                                         'mu_c': fixed(-1.0),\n",
    "                                         'sigma_g': fixed(0.5),\n",
    "                                         'sigma_c': fixed(value=0.5),\n",
    "                                         'plot_utility_row': fixed(False)})\n",
    "    display(widget_ui, widget_out)\n",
    "\n",
    "\n",
    "global global_plot_bayes_loss_utility_switcher\n",
    "global_plot_bayes_loss_utility_switcher = False\n",
    "def plot_bayes_loss_utility_switcher(what_to_plot):\n",
    "    global global_plot_bayes_loss_utility_switcher\n",
    "    if global_plot_bayes_loss_utility_switcher:\n",
    "        clear_output()\n",
    "    else:\n",
    "        global_plot_bayes_loss_utility_switcher = True\n",
    "\n",
    "    if what_to_plot == \"Gaussian\":\n",
    "        loss_f_options = Dropdown(\n",
    "                      options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "                      value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "        mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_prior\", continuous_update=True)\n",
    "        sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_prior\", continuous_update=True)\n",
    "        mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_true\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([loss_f_options, mu1_slider, sigma1_slider]),\n",
    "                          VBox([mu_true_slider, mu2_slider, sigma2_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_bayes_loss_utility_gaussian,\n",
    "                                                {'mu1': mu1_slider,\n",
    "                                                'mu2': mu2_slider,\n",
    "                                                'sigma1': sigma1_slider,\n",
    "                                                'sigma2': sigma2_slider,\n",
    "                                                'mu_true': mu_true_slider,\n",
    "                                                'loss_f': loss_f_options})\n",
    "\n",
    "    elif what_to_plot == \"Uniform\":\n",
    "        loss_f_options = Dropdown(\n",
    "                      options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "                      value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_true\", continuous_update=True)\n",
    "\n",
    "        widget_ui = HBox([VBox([loss_f_options, mu_slider, sigma_slider]),\n",
    "                          VBox([mu_true_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_bayes_loss_utility_uniform,\n",
    "                                        {'mu': mu_slider,\n",
    "                                        'sigma': sigma_slider,\n",
    "                                        'mu_true': mu_true_slider,\n",
    "                                        'loss_f': loss_f_options})\n",
    "\n",
    "    elif what_to_plot == \"Gamma\":\n",
    "\n",
    "        loss_f_options = Dropdown(\n",
    "                      options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "                      value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "\n",
    "        alpha_slider = FloatSlider(min=1.0, max=10.0, step=0.1, value=2.0, description=\"α_prior\", continuous_update=True)\n",
    "        beta_slider = FloatSlider(min=0.5, max=2.0, step=0.01, value=1.0, description=\"β_prior\", continuous_update=True)\n",
    "        # offset_slider = FloatSlider(min=-6.0, max=2.0, step=0.1, value=0.0, description=\"offset\", continuous_update=True)\n",
    "        offset_slider = fixed(0.0)\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_true\", continuous_update=True)\n",
    "        gaus_label = Label(value=\"normal likelihood\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "        gamma_label = Label(value=\"gamma prior\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "        widget_ui = HBox([VBox([loss_f_options, gamma_label, alpha_slider, beta_slider]),\n",
    "                          VBox([mu_true_slider, gaus_label, mu_slider, sigma_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_bayes_loss_utility_gamma,\n",
    "                                        {'alpha': alpha_slider,\n",
    "                                         'beta': beta_slider,\n",
    "                                         'offset': offset_slider,\n",
    "                                         'mu': mu_slider,\n",
    "                                         'sigma': sigma_slider,\n",
    "                                         'mu_true': mu_true_slider,\n",
    "                                         'loss_f': loss_f_options})\n",
    "\n",
    "    elif what_to_plot == \"Mixture of Gaussians\":\n",
    "        mu_m1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_mix_p\", continuous_update=True)\n",
    "        mu_m2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_mix_q\", continuous_update=True)\n",
    "        sigma_m1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_mix_p\", continuous_update=True)\n",
    "        sigma_m2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_mix_q\", continuous_update=True)\n",
    "        factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description=\"π\", continuous_update=True)\n",
    "        mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "        sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "        mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_true\", continuous_update=True)\n",
    "        loss_f_options = Dropdown(\n",
    "            options=[\"Mean Squared Error\", \"Absolute Error\", \"Zero-One Loss\"],\n",
    "            value=\"Mean Squared Error\", description=\"Loss: \")\n",
    "        empty_label = Label(value=\" \")\n",
    "\n",
    "        widget_ui = HBox([VBox([loss_f_options, mu_m1_slider, sigma_m1_slider]),\n",
    "                          VBox([mu_true_slider, mu_m2_slider, sigma_m2_slider]),\n",
    "                          VBox([empty_label, mu_slider, sigma_slider])])\n",
    "\n",
    "        widget_out = interactive_output(plot_bayes_loss_utility_mixture,\n",
    "                                        {'mu_m1': mu_m1_slider,\n",
    "                                         'mu_m2': mu_m2_slider,\n",
    "                                         'sigma_m1': sigma_m1_slider,\n",
    "                                         'sigma_m2': sigma_m2_slider,\n",
    "                                         'factor': factor_slider,\n",
    "                                         'mu': mu_slider,\n",
    "                                         'sigma': sigma_slider,\n",
    "                                         'mu_true': mu_true_slider,\n",
    "                                         'loss_f': loss_f_options})\n",
    "    display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title Helper Functions\n",
    "\n",
    "def gaussian(x, μ, σ):\n",
    "    \"\"\" Compute Gaussian probability density function for given value of the\n",
    "    random variable, mean, and standard deviation\n",
    "\n",
    "    Args:\n",
    "      x (scalar): value of random variable\n",
    "      μ (scalar): mean of Gaussian\n",
    "      σ (scalar): standard deviation of Gaussian\n",
    "\n",
    "    Returns:\n",
    "      scalar: value of probability density function\n",
    "    \"\"\"\n",
    "    return np.exp(-((x - μ) / σ)**2 / 2) / np.sqrt(2 * np.pi * σ**2)\n",
    "\n",
    "\n",
    "def gamma_pdf(x, α, β):\n",
    "    return gamma_distribution.pdf(x, a=α, scale=1/β)\n",
    "\n",
    "\n",
    "def mvn2d(x, y, mu1, mu2, sigma1, sigma2, cov12):\n",
    "    mvn = multivariate_normal([mu1, mu2], [[sigma1**2, cov12], [cov12, sigma2**2]])\n",
    "    return mvn.pdf(np.dstack((x, y)))\n",
    "\n",
    "\n",
    "def product_guassian(mu1, mu2, sigma1, sigma2):\n",
    "    J_1, J_2 = 1/sigma1**2, 1/sigma2**2\n",
    "    J_3 = J_1 + J_2\n",
    "    mu_prod = (J_1*mu1/J_3) + (J_2*mu2/J_3)\n",
    "    sigma_prod = np.sqrt(1/J_3)\n",
    "    return mu_prod, sigma_prod\n",
    "\n",
    "\n",
    "def reverse_product(mu3, sigma3, mu1, mu2):\n",
    "    J_3 = 1/sigma3**2\n",
    "    J_1 = J_3 * (mu3 - mu2) / (mu1 - mu2)\n",
    "    J_2 = J_3 * (mu3 - mu1) / (mu2 - mu1)\n",
    "    sigma1, sigma2 = 1/np.sqrt(J_1), 1/np.sqrt(J_2)\n",
    "    return sigma1, sigma2\n",
    "\n",
    "\n",
    "def calc_mean_mode_median(x, y):\n",
    "    \"\"\"\n",
    "\n",
    "    \"\"\"\n",
    "    pdf = y * (x[1] - x[0])\n",
    "    # Calc mode of an arbitrary function\n",
    "    mode = x[np.argmax(pdf)]\n",
    "\n",
    "    # Calc mean of an arbitrary function\n",
    "    mean = np.multiply(x, pdf).sum()\n",
    "\n",
    "    # Calc median of an arbitrary function\n",
    "    cdf = np.cumsum(pdf)\n",
    "    idx = np.argmin(np.abs(cdf - 0.5))\n",
    "    median = x[idx]\n",
    "\n",
    "    return mean, median, mode\n",
    "\n",
    "\n",
    "def calc_expected_loss(loss_f, posterior, x):\n",
    "    dx = x[1] - x[0]\n",
    "    expected_loss = np.zeros_like(x)\n",
    "    for i in np.arange(x.shape[0]):\n",
    "        loss = calc_loss_func(loss_f, x[i], x) # or mse or zero_one_loss\n",
    "        expected_loss[i] = np.sum(loss * posterior) * dx\n",
    "    return expected_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 1: Astrocat!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "<img src=\"https://github.com/NeuromatchAcademy/course-content/blob/master/tutorials/static/astrocat.png?raw=True\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Let's say you are a cat astronaut - Astrocat! You are navigating around space using a jetpack and your goal is to chase a mouse. \n",
    "\n",
    "Since you are a cat, you don't have opposable thumbs and cannot control your own jet pack. It can only be controlled by ground control back on Earth. \n",
    "\n",
    "For them to be able to guide you, they need to know where you are. They are trying to figure out your location. They have prior knowledge of your location - they know you like to hang out near the space mouse. They can also get an unreliable quick glimpse: they are taking a measurement of the hidden state of your location.\n",
    "\n",
    "They will try to figure out your location using Bayes rule and Bayesian decisions - as we will see throughout this tutorial.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "Astrocat is in space and we are considering the position along one dimension. So, the hidden state, s, is the true location. The satellites represent potential loss, and the space mouse, gain. Using indirect measurements, you as ground control, can estimate where Astrocat is or decide where it's likely better to send Astrocat.\n",
    "\n",
    "Remember, in this example, you can think of yourself as a scientist trying to decide where we belive Astrocat is, how to select a point estimate (single guess of location) based on possible errors, and how to account for the uncertainty we have about the location of the satellite and the space mouse. In fact, this is the kind of problem real scientists working to control remote satellites face! However, we can also think of this as what your brain has does when it wants to determine a target to make a movement or hit a tennis ball! A number of classic experiments use this kind of framing to study how *optimal* human decisions or movements are! Some examples are in the further reading document."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 2: Probability distribution of Astrocat location\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "We are going to think first about how Ground Control should estimate his position. We won't consider measurements yet, just how to represent the uncertainty we might have in general. We are now dealing with a continuous distribution - Astrocat's location can be any real number. In the last tutorial, we were dealing with a discrete distribution - the fish were either on one side or the other. \n",
    "\n",
    "So how do we represent the probability of each possible point (an infinite number) where the Astrocat could be? \n",
    "The Bayesian approach can be used on any probability distribution. While many variables in the world require representation using complex or unknown (e.g. empirical) distributions, we will be using the Gaussian distributions or extensions of it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 2.1: The Gaussian distribution\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "One distribution we will use throughout this tutorial is the **Gaussian distribution**, which is also sometimes called the normal distribution. \n",
    "\n",
    "This is a special, and commonly used, distribution for a couple reasons. It is actually the focus of one of the most important theorems in statistics: the Central Limit Theorem. This theorem tells us that if you sum a large number of samples of a variable, that sum is normally distributed *no matter what* the original distribution over a variable was. This is a bit too in-depth for us to get into now but check out links in the Bonus for more information. Additionally, Gaussians have some really nice mathematical properties that permit simple closed-form solutions to several important problems. As we will see later in this tutorial, we can extend Gaussians to be even more flexible and well approximate other distributions using mixtures of Gaussians. In short, the Gaussian is probably the most important continous distribution to understand and use.\n",
    "\n",
    "\n",
    "Gaussians have two parameters. The **mean** $\\mu$, which sets the location of its center. Its \"scale\" or spread is controlled by its **standard deviation** $\\sigma$ or its square, the **variance** $\\sigma^2$. These can be a bit easy to mix up: make sure you are careful about whether you are referring to/using standard deviation or variance.\n",
    "\n",
    "The equation for a Gaussian distribution on a variable $x$ is:\n",
    "\n",
    "$$\n",
    "\\mathcal{N}(\\mu,\\sigma^2) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp\\left(\\frac{-(x-\\mu)^2}{2\\sigma^2}\\right)\n",
    "$$\n",
    "\n",
    "In our example, $x$ is the location of the Astrocat in one direction. $\\mathcal{N}(\\mu,\\sigma^2)$ is a standard notation to refer to a **N**ormal (Gaussian) distribution. For example, $\\mathcal{N}(0, 1)$ denotes a Gaussian distribution with mean 0 and variance 1. The exact form of this equation is not particularly intuitive, but we will see how mean and standard deviation values affect the probability distribution.\n",
    "\n",
    "\n",
    "We won't implement a Gaussian distribution in code here but please refer to the pre-reqs refresher W0D5 T1 to do this if you need further clarification.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the function `gaussian`\n",
    "def gaussian(x, μ, σ):\n",
    "    return np.exp(-((x - μ) / σ)**2 / 2) / np.sqrt(2 * np.pi * σ**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 2.1: Exploring Gaussian parameters:\n",
    "\n",
    "Let's explore how the parameters of a Gaussian affect the distribution. Play with the demo below by changing the mean $\\mu$ and standard deviation $\\sigma$.\n",
    "\n",
    "Discuss the following:\n",
    "\n",
    "1. What does increasing $\\mu$ do? What does increasing $\\sigma$ do?\n",
    "2. If you wanted to know the probability of an event happing at $0$, can you find two different $\\mu$ and $\\sigma$ values that produce the same probabilty of an event at $0$?\n",
    "3. How many Gaussian's could produce the same probabilty at $0$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "widget = interact(plot_gaussian,\n",
    "                     μ = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0, continuous_update=False),\n",
    "                     σ = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, continuous_update=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_66e77019.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 2.2: Multiplying Gaussians\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "When we multiply Gaussians, we are not multiplying random variables but the actual underlying distributions. If we multiply two Gaussian distributions, with means $\\mu_1$ and $\\mu_2$ and standard deviations $\\sigma_1$ and $\\sigma_2$, we get another Gaussian. The Gaussian resulting from the multiplication will have mean $\\mu_3$ and standard deviation $\\sigma_3$ where:\n",
    "\n",
    "$$\n",
    "\\mu_{3} = a\\mu_{1} + (1-a)\\mu_{2}\n",
    "$$\n",
    "$$\n",
    "\\sigma_{3}^{-2} = \\sigma_{1}^{-2} + \\sigma_{2}^{-2}\\\\\n",
    "a = \\frac{\\sigma_{1}^{-2}}{\\sigma_{1}^{-2} + \\sigma_{2}^{-2}}\n",
    "$$\n",
    "\n",
    "This may look confusing but keep in mind that the information in a Gaussian is the inverse of its variance: $\\frac{1}{\\sigma^2}$. Basically, when multiplying Gaussians, the mean of the resulting Gaussian is a weighted average of the original means, where the weights are proportional to the amount of information of that Gaussian. The information in the resulting Gaussian is equal to the sum of informations of the original two. We'll dive into this in the next demo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 2.2: Multiplying Gaussians \n",
    "\n",
    "We have implemented the multiplication of two Gaussians for you. Using the following widget, we are going to think about the information and combination of two Gaussians. It In our case, imagine we want to find the middle location between the satellite and the space mouse. This would be the center (average) of the two locations. Because we have uncertainty, we need to weight our uncertainty in thinking about the most likely place. \n",
    "\n",
    "In this demo, $\\mu_{1}$ and $\\sigma_{1}$ are the mean and standard deviation of the distribution over satellite location, $\\mu_{2}$ and $\\sigma_{2}$ are the mean and standard deviation of the distribution over space mouse location, and $\\mu_{3}$ and $\\sigma_{3}$ are the mean and standard deviation of the distribution over the center location (gained by multiplying the first two).\n",
    "\n",
    "Questions:\n",
    "\n",
    "1. What is your uncertainty (how much information) do you have about $\\mu_{3}$ with $\\mu_{1} = -2, \\mu_{2} = 2, \\sigma_{1} = \\sigma_{2} = 0.5$?\n",
    "2. What happens to your estimate of $\\mu_{3}$ as $\\sigma_{2} \\to \\infty$? (In this case, $\\sigma$ only goes to 11... but that should be loud enough.)\n",
    "3. What is the difference in your estimate of $\\mu_{3}$ if $\\sigma_{1} = \\sigma_{2} = 11$? What has changed from the first example?\n",
    "4. Set $\\mu_{1} = -4, \\mu_{2} = 4$ and change the $\\sigma$s so that $\\mu_{3}$ is close to $2$. How many $\\sigma$s will produce the same $\\mu_{3}$?\n",
    "5. Continuing, if you set $\\mu_{1} = 0$, what $\\sigma$ do you need to change so $\\mu_{3} ~= 2$?\n",
    "6. If $\\sigma_{1} = \\sigma_{2} = 0.1$, how much information do you have about the average?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "mu1_slider = FloatSlider(min=-5.0, max=-0.51, step=0.01, value=-2.0, description=\"µ_1\",continuous_update=True)\n",
    "mu2_slider = FloatSlider(min=0.5, max=5.01, step=0.01, value=2.0, description=\"µ_2\",continuous_update=True)\n",
    "sigma1_slider = FloatSlider(min=0.1, max=11.01, step=0.01, value=1.0, description=\"σ_1\", continuous_update=True)\n",
    "sigma2_slider = FloatSlider(min=0.1, max=11.01, step=0.01, value=1.0, description=\"σ_2\", continuous_update=True)\n",
    "distro_1_label = Label(value=\"Satellite\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "distro_2_label = Label(value=\"Space Mouse\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "widget_ui = HBox([VBox([distro_1_label, mu1_slider, sigma1_slider]),\n",
    "                  VBox([distro_2_label, mu2_slider, sigma2_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_information, {'mu1': mu1_slider,\n",
    "                                                    'mu2': mu2_slider,\n",
    "                                                    'sigma1': sigma1_slider,\n",
    "                                                    'sigma2': sigma2_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_9b8c9143.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "To start thinking about how we might use these concepts directly in systems neuroscience, imagine you want to know how much information is gained combining (averaging) the response of two neurons that represent locations in sensory space (think: how much information is shared by their receptive fields). You would be multiplying Gaussians!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 2.3: Mixtures of Gaussians\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "What if our continuous distribution isn't well described by a single bump? For example, what if the Astrocat is often either in one place or another - a Gaussian distribution would not capture this well! We need a multimodal distribution. Luckily, we can extend Gaussians into a *mixture of Gaussians*, which are more complex distributions.  \n",
    "\n",
    "In a Gaussian mixture distribution, you are essentially adding two or more weighted standard Gaussian distributions (and then normalizing so everything integrates to 1). Each standard Gaussian involved is described, as normal, by its mean and standard deviation. Additional parameters in a mixture of Gaussians are the weights you put on each Gaussian (π). The following demo should help clarify how a mixture of Gaussians relates to the standard Gaussian components. We will not cover the derivation here but you can work it out as a bonus exercise.\n",
    "\n",
    "Mixture distributions are a common tool in Bayesian modeling and an important tool in general.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 2.3: Exploring Gaussian mixtures\n",
    "\n",
    "We will examing a mixture of two Gaussians. We will have one weighting parameter, π, that tells us how to weight one of the Gaussians. The other is weighted by 1 - π. \n",
    "\n",
    "Use the following widget to experiment with the parameters of each Gaussian and the mixing weight ($\\pi$) to undersand how the mixture of Gaussian distribution behaves.\n",
    "\n",
    "Discuss the following questions:\n",
    "\n",
    "1. What does increasing the weight $\\pi$ do to the mixture distribution (dark blue)? \n",
    "2. How can you make the two bumps of the mixture distribution further apart?\n",
    "3. Can you make the mixture distribution have just one bump (like a Gaussian)?\n",
    "4. Any other shapes you can make the mixture distribution resemble other than one nicely rounded bump or two separate bumps?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-1.0, description=\"µ_p\", continuous_update=True)\n",
    "mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=1.0, description=\"µ_q\", continuous_update=True)\n",
    "sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_p\", continuous_update=True)\n",
    "sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_q\", continuous_update=True)\n",
    "factor_slider = FloatSlider(min=0.1, max=0.9, step=0.01, value=0.5, description=\"π\", continuous_update=True)\n",
    "distro_1_label = Label(value=\"p(x)\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "distro_2_label = Label(value=\"q(x)\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "mixture_label = Label(value=\"mixing coefficient\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "widget_ui = HBox([VBox([distro_1_label, mu1_slider, sigma1_slider]),\n",
    "                  VBox([distro_2_label, mu2_slider, sigma2_slider]),\n",
    "                  VBox([mixture_label, factor_slider])])\n",
    "\n",
    "widget_out = interactive_output(gaussian_mixture, {'mu1': mu1_slider,\n",
    "                                                    'mu2': mu2_slider,\n",
    "                                                    'sigma1': sigma1_slider,\n",
    "                                                    'sigma2': sigma2_slider,\n",
    "                                                    'factor': factor_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_93a3e0b7.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 3: Utility \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "We want to know where Astrocat is. If we were asked to provide the coordinates, for example to display them for Ground Control or to note them in a log, we are not going to provide the whole probability distribution! We will give a single set of coordinates, but we first need to estimate those coordinates. Just like in the last tutorial, this may not be as easy as just what is most likely: we want to know how good or bad it is if we guess a certain location and the Astrocat is in another. \n",
    "\n",
    "As we have seen, utility represents the gain (or if negative, loss) for us if we take a certain action for a certain value of the hidden state. In our continuous example, we need a function to be able to define the utility with respect to all possible continuous values of the state. Our action here is our guess of the Astrocat location.\n",
    "\n",
    "We are going to explore this for the Gaussian distribution, where our estimate is $\\hat{\\mu}$ and the true hidden state we are interested in is $\\mu$. \n",
    "\n",
    "A loss function determines the \"cost\" (or penalty) of estimating $\\hat \\mu$ when the true or correct quantity is really $\\mu$ (this is essentially the cost of the error between the true hidden state we are interested in: $\\mu$ and our estimate: $\\hat \\mu$). A loss function is equivalent to a negative utility function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 3.1: Standard loss functions\n",
    "\n",
    "\n",
    "There are lots of different possible loss functions. We will focus on three: **mean-squared error** where the loss is the different between truth and estimate squared, **absolute error** where the loss is the absolute difference between truth and estimate, and **Zero-one Loss** where the loss is 1 unless we're exactly right (the estimate equals the truth). We can represent these with the following formulas:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\textrm{Mean Squared Error} &=& (\\mu - \\hat{\\mu})^2 \\\\ \n",
    "\\textrm{Absolute Error} &=& \\big|\\mu - \\hat{\\mu}\\big| \\\\ \n",
    "\\textrm{Zero-One Loss} &=& \\begin{cases}\n",
    "                            0,& \\textrm{if } \\mu = \\hat{\\mu} \\\\\n",
    "                            1,              & \\textrm{otherwise}\n",
    "                            \\end{cases}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "We will now explore how these different loss functions change our expected utility!\n",
    "\n",
    "Check out the next cell to see the implementation of each loss in the function `calc_loss_func`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the function `calc_loss_func`\n",
    "\n",
    "def calc_loss_func(loss_f, mu_true, x):\n",
    "    error = x - mu_true\n",
    "    if loss_f == \"Mean Squared Error\":\n",
    "        loss = (error)**2\n",
    "    elif loss_f == \"Absolute Error\":\n",
    "        loss = np.abs(error)\n",
    "    elif loss_f == \"Zero-One Loss\":\n",
    "        loss = (np.abs(error) >= 0.03).astype(np.float)\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive demo 3: Exploring Loss with different distributions\n",
    "\n",
    "Let's see how our loss functions interact with probability distributions to affect expected utility and consequently, the action we take.\n",
    "\n",
    "Play with the widget below and discuss the following:\n",
    "\n",
    "1. With a Gaussian distribution, does the peak of the expected utility ever change position on the x-axis for the three different losses? This peak denotes the action we would choose (the location we would guess) so in other words, would the different choices of loss function affect our action?\n",
    "2. With a mixture of Gaussian distribution with two bumps, does the peak of the expected loss ever change position on the x-axis for the three different losses?\n",
    "3.  Find parameters for a mixture of Gaussians that results in the mean, mode, and median all being distinct (not equal to one another). With this distribution, how does the peak of the expected utility correspond to the mean/median/mode of the probability distribution for each of the three losses?\n",
    "4. When the mixture of Gaussians has two peaks that are exactly the same height, how many modes are there?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "widget = interact(loss_plot_switcher,\n",
    "                  what_to_plot = Dropdown(\n",
    "                      options=[\"Gaussian\", \"Mixture of Gaussians\"],\n",
    "                      value=\"Gaussian\", description=\"Distribution: \"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_79bc9081.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "You can see that what coordinates you would provide for Astrocat aren't necessarily easy to guess just from the probability distribution. You need the concept of utility/loss and a specific loss function to determine what estimate you should give.\n",
    "\n",
    "For symetric distributions, you will find that the mean, median and mode are the same. However, for distributions with *skew*, like the Gamma distribution or the Exponential distribution, these will be different. You will be able to explore more distributions as priors below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 3.2: A more complex loss function\n",
    "\n",
    "\n",
    "The loss functions we just explored were fairly simple and are often used. However, life can be complicated and in this case, Astrocat cares about both being near the space mouse and avoiding the satellites. This means we need a more complex loss function that captures this! \n",
    "\n",
    "We know that we want to estimate Astrocat to be closer to the mouse, which is safe and desirable, but further away from the satellites, which is dangerous! So, rather than thinking about the *Loss* function, we will consider a generalized utility function that considers gains and losses that *matter* to Astrocat!\n",
    "\n",
    "In this case, we can assume that depending on our uncertainty about Astrocat's probable location, we may want to 'guess' that Astrocat is close to 'good' parts of space and further from 'bad' parts of space. We will model these utilities as Gaussian gain and loss regions--and we can assume the width of the Gaussian comes from our uncertainty over where the Space Mouse and satellite are.\n",
    "\n",
    "Let's explore how this works in the next interactive demo.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive demo 3.2: Complicated cat costs\n",
    "\n",
    "Now that we have explored *Loss* functions that can be used to determine both formal *estimators* and our expected loss given our error, we are going to see what happens to our estimates if we use a generalized utility function.\n",
    "\n",
    "Questions:\n",
    "\n",
    "1. As you change the $\\mu$ of Astrocat, what happens to the expected utility?\n",
    "2. Can the EU be exactly zero everywhere?\n",
    "3. Can the EU be zero in a region around Astrocat but positive and negative elsewhere?\n",
    "4. As our uncertainty about Astrocat's position increases, what happens to the expected utility?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0, description=\"µ\", continuous_update=True)\n",
    "mu_g_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-1.0, description=\"µ_gain\", continuous_update=True)\n",
    "mu_c_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=1.0, description=\"µ_cost\", continuous_update=True)\n",
    "sigma_slider =  FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ\", continuous_update=True)\n",
    "sigma_g_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_gain\", continuous_update=True)\n",
    "sigma_c_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_cost\", continuous_update=True)\n",
    "\n",
    "distro_label = Label(value=\"probability\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "gain_label = Label(value=\"gain\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "loss_label = Label(value=\"loss\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "\n",
    "widget_ui = HBox([VBox([distro_label, mu_slider, sigma_slider]),\n",
    "                  VBox([gain_label, mu_g_slider, sigma_g_slider]),\n",
    "                  VBox([loss_label, mu_c_slider, sigma_c_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_simple_utility_gaussian,\n",
    "                                {'mu': mu_slider,\n",
    "                                'sigma': sigma_slider,\n",
    "                                'mu_g': mu_g_slider,\n",
    "                                'mu_c': mu_c_slider,\n",
    "                                'sigma_g': sigma_g_slider,\n",
    "                                 'sigma_c': sigma_c_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_96977d91.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 4: Correlation and marginalization\n",
    "\n",
    "\n",
    "*Estimated timing to here from start of tutorial: 40 min*\n",
    "\n",
    "In this section we will explore a two dimensional Gaussian, often defined as a two-dimension vector of Gaussian random variables. This is, in essence, the joint distribution of two Gaussian random variables.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 4.1: Correlation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "If the two variables in a two dimensional Gaussian are independent, looking at one tells us nothing about the other. But what if the the two variables are correlated (covary)?\n",
    "\n",
    "The covariance of two Gaussians with means $\\mu_X$ and $\\mu_Y$ and standard deviations $\\sigma_X$ and $\\sigma_Y$is:\n",
    "\n",
    "$$\n",
    "\\sigma_{XY} = E[(X-\\mu_{X})(Y-\\mu_{Y})]\n",
    "$$\n",
    "\n",
    "$E$ here denotes the expected value. So the covariance is the expected value of the random variable X minus the mean of the Gaussian distribution on X times the random variable Y minus the mean of the Gaussian distribution on Y.\n",
    "\n",
    "The correlation is the covariance normalized, so that it goes between -1 (exactly anticorrelated) to 1 (exactly correlated).\n",
    "\n",
    "$$\n",
    "\\rho_{XY} = \\frac{\\sigma_{XY}}{\\sigma_{X}\\sigma_{Y}}\n",
    "$$\n",
    "\n",
    "These are key concepts and while we are considering two hidden states (or two random variables), they extend to $N$ dimensional vectors of Gaussian random variables. You will find these used all over computational neuroscience.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive demo 4.1: Covarying 2D Gaussian\n",
    "\n",
    "Let's explore this 2D Gaussian (i.e. joint distribution of two Gaussians).\n",
    "\n",
    "Use the following widget to think about the following questions:\n",
    "\n",
    "1. If these variables represent hidden states we care about, what does observing one tell us about the other? How does this depend on the correlation?\n",
    "2. How does the shape of the distribution change when we change the means? The variances? The correlation?\n",
    "3. If we want to isolate one or the other hidden state distributions, what do we need to do? (Hint: think about Tutorial 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute the cell to enable the widget\n",
    "\n",
    "mu1_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"µ_x\", continuous_update=True)\n",
    "# mu2_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"µ_y\", continuous_update=True, orientation='vertical')\n",
    "mu2_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"µ_y\", continuous_update=True)\n",
    "sigma1_slider = FloatSlider(min=0.1, max=1.5, step=0.01, value=0.5, description=\"σ_x\", continuous_update=True)\n",
    "# sigma2_slider = FloatSlider(min=0.1, max=1.5, step=0.01, value=0.5, description=\"σ_y\", continuous_update=True, orientation='vertical')\n",
    "sigma2_slider = FloatSlider(min=0.1, max=1.5, step=0.01, value=0.5, description=\"σ_y\", continuous_update=True)\n",
    "corr_slider = FloatSlider(min=-0.99, max=0.99, step=0.01, value=0.0, description=\"ρ\", continuous_update=True)\n",
    "\n",
    "distro1_label = Label(value=\"p(x)\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "distro2_label = Label(value=\"p(y)\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "corr_label = Label(value=\"correlation\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "\n",
    "\n",
    "widget_ui = HBox([VBox([distro1_label, mu1_slider, sigma1_slider]),\n",
    "                  VBox([corr_label, corr_slider]),\n",
    "                  VBox([distro2_label, mu2_slider, sigma2_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_mvn2d, {'mu1': mu1_slider,\n",
    "                                            'mu2': mu2_slider,\n",
    "                                            'sigma1': sigma1_slider,\n",
    "                                            'sigma2': sigma2_slider,\n",
    "                                            'corr': corr_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_159a3f69.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 4.2: Marginalization and information\n",
    "\n",
    "\n",
    "We learned in Tutorial 1 that if we want to measure the probability of one or another variable, we need to average over the other. When we extend this to the correlated Gaussians we just played with, marginalization works the same way. Let's say that the two variables reflect Astrocat's position in space (in two dimensions). If we want to get our uncertainty about Astrocat's X or Y position, we need to marginalize. \n",
    "\n",
    "However, let's imagine we have a measurement from one of the variables, for example X, and we want to understand the uncertainty we have in Y. We no longer want to marginalize because we know X, we don't need to ignore it! Instead, we can calculate the conditional probability $P(Y|X=x)$. You will explore the relationship between these two concepts in the following interactive demo.\n",
    "\n",
    "But first, let's remember that we can also think about the amount of uncertainty as inversely proportional to the amount of information we have about each variable. This is important, because the joint information is determined by the correlation. For our Bayesian approach, the important intuition is that we can also think about the mutual information between the prior and the likelihood following a measurement.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive demo 4.2: Marginalizing 2D Gaussians\n",
    "\n",
    "Use the following widget to think consider the following questions:\n",
    "\n",
    "1. When is the marginal distribution the same as the conditional probability distribution? Why?\n",
    "2. If $\\rho$ is large, how much information can we gain (in addition) looking at both variables vs just considering one?\n",
    "3. If $\\rho$ is close to zero, but the variances of the two variables are very different, what happens to the conditional probability compared to the marginals? As $\\rho$ changes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "\n",
    "sigma1_slider = FloatSlider(min=0.1, max=1.1, step=0.01, value=0.5, description=\"σ_x\", continuous_update=True)\n",
    "sigma2_slider = FloatSlider(min=0.1, max=1.1, step=0.01, value=0.5, description=\"σ_y\", continuous_update=True)\n",
    "c_x_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"Cx\", continuous_update=True)\n",
    "c_y_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"Cy\", continuous_update=True)\n",
    "corr_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0, description=\"ρ\", continuous_update=True)\n",
    "\n",
    "distro1_label = Label(value=\"x\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "distro2_label = Label(value=\"y\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "corr_label = Label(value=\"correlation\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "\n",
    "widget_ui = HBox([VBox([distro1_label, sigma1_slider, c_x_slider]),\n",
    "                  VBox([corr_label, corr_slider]),\n",
    "                  VBox([distro2_label, sigma2_slider, c_y_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_marginal, {'sigma1': sigma1_slider,\n",
    "                                                'sigma2': sigma2_slider,\n",
    "                                                'c_x': c_x_slider,\n",
    "                                                'c_y': c_y_slider,\n",
    "                                                'corr': corr_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_3ec1ad5b.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 5: Bayes' theorem for continuous distributions\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "The continous case allows us to consider how the shape of the posterior distribution can differ from the prior. The Gaussian case is the most fundemental, but asymetric priors (or likelihoods) and posteriors allow us to see how the mean, median and mode can be effected differently when we apply Bayes' theorem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 5.1: The Gaussian example\n",
    "\n",
    "Bayes' rule tells us how to combine two sources of information: the prior (e.g., a noisy representation of Ground Control's expectations about where Astrocat is) and the likelihood (e.g., a noisy representation of the Astrocat after taking a measurement), to obtain a posterior distribution (our belief distribution) taking into account both pieces of information. Remember Bayes' rule:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\text{Posterior} = \\frac{ \\text{Likelihood} \\times \\text{Prior}}{ \\text{Normalization constant}}\n",
    "\\end{eqnarray}\n",
    "\n",
    "We will look at what happens when both the prior and likelihood are Gaussians. In these equations, $\\mathcal{N}(\\mu,\\sigma^2)$ denotes a Gaussian distribution with parameters $\\mu$ and $\\sigma^2$:\n",
    "$$\n",
    "\\mathcal{N}(\\mu, \\sigma) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} \\; \\exp \\bigg( \\frac{-(x-\\mu)^2}{2\\sigma^2} \\bigg)\n",
    "$$\n",
    "\n",
    "\n",
    "When both the prior and likelihood are Gaussians, Bayes Rule translates into the following form:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "\\text{Likelihood} &=& \\mathcal{N}(\\mu_{likelihood},\\sigma_{likelihood}^2) \\\\\n",
    "\\text{Prior} &=& \\mathcal{N}(\\mu_{prior},\\sigma_{prior}^2) \\\\\n",
    "\\text{Posterior} &=& \\mathcal{N}\\left( \\frac{\\sigma^2_{likelihood}\\mu_{prior}+\\sigma^2_{prior}\\mu_{likelihood}}{\\sigma^2_{likelihood}+\\sigma^2_{prior}}, \\frac{\\sigma^2_{likelihood}\\sigma^2_{prior}}{\\sigma^2_{likelihood}+\\sigma^2_{prior}} \\right) \\\\\n",
    "&\\propto& \\mathcal{N}(\\mu_{likelihood},\\sigma_{likelihood}^2) \\times \\mathcal{N}(\\mu_{prior},\\sigma_{prior}^2)\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "We get the parameters of the posterior from multiplying the Gaussians, just as we did in Secton 2.2.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 5.1: Gaussian Bayes\n",
    "Let's consider the following questions using the following interactive demo:\n",
    "\n",
    "1. For a Gaussian posterior, explain how the information seems to be combining. (Hint: think about the prior exercises!)\n",
    "2. What is the difference between the posterior here and the Gaussian that represented the average of two Gaussians in the exercise above?\n",
    "3. How should we think about the relative weighting of information between the prior and posterior?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_prior\", continuous_update=True)\n",
    "mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_prior\", continuous_update=True)\n",
    "sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "distro1_label = Label(value=\"prior distribution\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "distro2_label = Label(value=\"likelihood distribution\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "\n",
    "widget_ui = HBox([VBox([distro1_label, mu1_slider, sigma1_slider]),\n",
    "                  VBox([distro2_label, mu2_slider, sigma2_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_bayes, {'mu1': mu1_slider,\n",
    "                                            'mu2': mu2_slider,\n",
    "                                            'sigma1': sigma1_slider,\n",
    "                                            'sigma2': sigma2_slider})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_db1f6cb5.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 5.2: Exploring priors\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 5.2: Prior exploration\n",
    "\n",
    "What would happen if we had a different prior distribution for Astrocat's location? Bayes' Rule works exactly the same way if our prior is not a Guassian (though the analytical solution may be far more complex or impossible). Let's look at how the posterior behaves if we have a different prior over Astrocat's location.\n",
    "\n",
    "Consider the following questions:\n",
    "\n",
    "1. Why does the posterior not look Gaussian when you use a non-Gaussian prior?\n",
    "2. What does having a flat prior mean?\n",
    "3. How does the Gamma prior behave differently than the others?\n",
    "4. From what you know, can you imagine the likelihood being something other than a Gaussian?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "widget = interact(plot_prior_switcher,\n",
    "                  what_to_plot = Dropdown(\n",
    "                      options=[\"Gaussian\", \"Mixture of Gaussians\", \"Uniform\", \"Gamma\"],\n",
    "                      value=\"Gaussian\", description=\"Prior: \"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_9753d295.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 6: Bayesian decisions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "Bayesian decisions in continuous dimensions are the same as we saw in Tutorial 1 for the binary case. The only difference is that now, our Expected Utility is calculated using an integral and all of our probability distributions are continuous."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "## Section 6.1: Bayesian estimation on the posterior\n",
    "\n",
    "Now that we understand that the posterior can be something other than a Gaussian, let's revisit **Loss** functions. In this case, we can see that the posterior can take many forms. \n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 6.1: Standard loss functions with various priors\n",
    "\n",
    "Questions:\n",
    "\n",
    "1. If we have a bi-modal prior, how do the different loss functions potentially inform us differently about what we learn?\n",
    "2. Why do the different loss functions behavior differently with respect to the shape of the posterior? When do they produce different expected loss?\n",
    "3. For the mixture of Gaussians, describe the situations where the expected loss will look different from the Gaussian case.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "widget = interact(plot_bayes_loss_utility_switcher,\n",
    "                  what_to_plot = Dropdown(\n",
    "                      options=[\"Gaussian\", \"Mixture of Gaussians\", \"Uniform\", \"Gamma\"],\n",
    "                      value=\"Gaussian\", description=\"Prior: \"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_d9631fab.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Section 6.2: Bayesian decisions\n",
    "\n",
    "Finally, we can combine everything we have learned so far! \n",
    "\n",
    "Now, let's imagine we have just received a new measurement of Astrocat's location. We need to think about how we want to decide where Astrocat is, so that we can decide how far to tell Astrocat to move. However, we want to account for the satellite and Space Mouse location in this estimation. If we make an error towards the satellite, it's worse than towards Space Mouse. So, we will use our more complex utility function from Section 3.2. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Interactive Demo 6.2: Complicated cat costs with various priors\n",
    "\n",
    "\n",
    "Questions:\n",
    "\n",
    "1. If you have a weak prior and likelihood, how much are you relying on the utility function to guide your estimation?\n",
    "2. If you get a good measurement, that is a likelihood with low variance, how much does this help?\n",
    "3. Which of the factors are most important in making your decision?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description=\"µ_prior\", continuous_update=True)\n",
    "mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description=\"µ_likelihood\", continuous_update=True)\n",
    "sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_prior\", continuous_update=True)\n",
    "sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_likelihood\", continuous_update=True)\n",
    "mu_g_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0, description=\"µ_gain\", continuous_update=True)\n",
    "sigma_g_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description=\"σ_gain\", continuous_update=True)\n",
    "\n",
    "dist_label = Label(value=\"loss distance: µ1_c - µ2_c\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "loc_label = Label(value=\"loss center: (µ1_c + µ2_c) / 2\", layout=Layout(display=\"flex\", justify_content=\"center\"))\n",
    "mu_dist_slider = FloatSlider(min=0.0, max=8.0, step=0.01, value=4.0, description=\"distance\", continuous_update=True)\n",
    "mu_loc_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0, description=\"center\", continuous_update=True)\n",
    "\n",
    "widget_ui = HBox([VBox([mu1_slider, sigma1_slider, mu2_slider, sigma2_slider]),\n",
    "                  VBox([dist_label, mu_dist_slider, loc_label, mu_loc_slider]),\n",
    "                  VBox([mu_g_slider, sigma_g_slider])])\n",
    "\n",
    "widget_out = interactive_output(plot_utility_mixture_dist,\n",
    "                                    {'mu1': mu1_slider,\n",
    "                                    'sigma1': sigma1_slider,\n",
    "                                    'mu2': mu2_slider,\n",
    "                                    'sigma2': sigma2_slider,\n",
    "                                    'mu_g': mu_g_slider,\n",
    "                                    'sigma_g': sigma_g_slider,\n",
    "                                    'mu_dist': mu_dist_slider,\n",
    "                                    'mu_loc': mu_loc_slider,\n",
    "                                    'plot_utility_row': fixed(True)})\n",
    "display(widget_ui, widget_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W7_BayesianDecisions/solutions/W7_Tutorial2_Solution_9b83b02b.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "In this tutorial, you extended your exploration of Bayes Rule and the Bayesian approach in the context of finding and choosling a location for Astrocat.\n",
    "\n",
    "Specifically, we covered:\n",
    "\n",
    "* The Gaussian distribution and its properties\n",
    "\n",
    "* That the likelihood is the probability of the measurement given some hidden state\n",
    "\n",
    "* Information shared between Gaussians (via multiplication of PDFs and via two-dimensional distributions)\n",
    "\n",
    "* That how the prior and likelihood interact to create the posterior, the probability of the hidden state given a measurement, depends on how they covary\n",
    "\n",
    "* That utility is the gain from each action and state pair, and the expected utility for an action is the sum of the utility for all state pairs, weighted by the probability of that state happening. You can then choose the action with highest expected utility.\n",
    " "
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W7_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "interpreter": {
   "hash": "9516f62da91337f10c2adbe814d9c63a4b08f8271333386358218606edb781e3"
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3.7.11 64-bit ('pw3': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
